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Abstract 

We study a model of flocking for a very large system (N=320,000) numeri- 
cally. We find that in the long wavelength, long time limit, the fluctuations of 
the velocity and density fields are carried by propagating sound modes, whose 
dispersion and damping agree quantitatively with the predictions of our pre- 
vious work using a continuum equation. We find that the sound velocity is 
anisotropic and characterized by its speed c for propagation perpendicular to 
the mean velocity < v >, < v > itself, and a third velocity A < v >, arising 
explicitly from the lack of Galilean invariance in flocks. 
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The dynamics of "flocking" behavior of living things, such as birds, fish, wildebeest, 
slime molds and bacteria has long attracted a great deal of attention among biologists, 
computer animators and physicists It is crucial to correctly describe the interaction 

between members of the flock in order to understand and model the flocking behavior. 
As summarized in [2], a large flock does not have a global leader; instead, the impressive 
collective flocking phenomena is caused by individual members of the flock following the 
motion of their neighbors. 

In our earlier work 0], we studied the flocking dynamics by using continuum equations 
for the coarse-grained density field p(x,t) and velocity field v(x,t), written as: 

d t v + X(v ■ V)v = a v - /3\v\ 2 v - VP + £> L V(V • v) 

+ D 1 V 2 v + D 2 (v- Vfv + f (1) 
H + V»)=0 (2) 

where /3, Di, D 2 and D L are all positive, and a < in the disordered phase and a > in 
the ordered state. The a and (3 terms simply make the local v have a non-zero magnitude 
(= Ja/0) in the ordered phase. P>l,\,i are diffusion constants. The Gaussian random noise 
/ has correlations: < fi{r 1 t)fj{r' 1 t l ) >= A5ij5 d (f — r')5(t — t') where A is a constant, and 
i , j denote Cartesian components. Finally, the pressure P = P(p) = X^i°Vt(p — Po) n , 
where po is the mean of the local number density p(r) and a n are coefficients in the pressure 
expansion. The final equation (2) reflects conservation of birds. 

In [4], we considered the special case of (1) with A = 1. Just as the absence of the 
Galilean invariance for the flock motion allows a and (3 ^ in equation (1), likewise A need 
not be = 1. In [5] and this paper, we consider the more generic case A ^ 1, which leads to 
a different direction dependence of the sound speed than when A = 1 |J. 

In the ordered phase where a > 0, the velocity field and the density field can be written 
as: v = v 3 x\\ + 5v , p = po + 6p where po and v s x\\ are the space averaged density and 
spontaneous symmetry broken velocity respectively. The spontaneous symmetry breaking 
of a vector field leads to large "Goldstone mode" fluctuations; in flocks, this mode is v±, 
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the projection of Sv perpendicular to x\\ (we will hereafter use "||" (" _L ") to denote the 
projection of any vector along (perpendicular to) afn). Indeed, for equilibrium systems, such 
fluctuations are strong enough in two dimensions to destroy the long range order J7[. One 
of the remarkable predictions of our continuum model of flocking is that the ordered state 
is stable even in two dimensions due to the non-equilibrium effect of the nonlinear terms. 
The mean squared fluctuations in Fourier space in 2D are: 

< > = fc^LtM (4) 

S(q,u) 

where the denominator S(q,u) = [(u — v s q\\)(w — ^v s q\\) — c 2 ^] 2 + [(u — ^s<?||)(-D;i(<?)<zi + 
-D||gjj) — \v s D p q^} 2 . c = v /cr 1 p , D p = c 2 /a, D\\ = D 1 + D 2 + D p and D±, the renormalized 
diffusion constant, scales as: 

Df(q ± , qil ;X,p ,a n )=qr 2 f(^) , (5) 

where the exponents z and ( are found by RG analysis to be z = | and ( = | for two 
dimensions, and the scaling function f(x) is universal up to an overall, nonuniversal q and 
u) independent scale factor. 

From the above expressions, the correlation functions will have peaks around ooo(q) which 
satisfies: (uj — v s q\\)(ujQ — \v s q\\) — c 2 q\ = with solutions: u = fl±(q) = |(1 + \)v s q\\ ± 
(|(1 — X) 2 v 2 q 2 + (^q 2 ^) 1 / 2 . This implies that for the wavevector (q\\,q±) = q(cos(9 q ),sm(6 q )) 
at small q , there should be two peaks in the power spectrum located around u = c±(9 q )q 
with the sound speeds: 

c±{6 q ) = ~(l + AKcos(0 g ) 

±(l(l-A) 2 ^co S 2 ^)+c 2 S in 2 (^))^ ( 6 ) 

The relative strength of the two peaks varies with 6 q . It is not hard to see from eqs. (3), 
(4), that at 9 q ~ , < \Sp(q, uo)\ 2 > will only have a peak with corresponding wave velocity 
c+(9 q = 0) = v s , and < \v±(q, u>)\ 2 > will only have a peak at C-(9 q = 0) = \v s . 



In this paper, we study a discrete model numerically to test the predictions made by our 
continuum theory. The model we use is very similar to the one studied by Vicsek et al j3]. 
Following [1], we call our simulated flockers "boids". At a given time t, the position and 
the direction of the velocity for each boid are given as (fi(t) ,9i(t)) for i = 1, 2, ■ • • , N. The 
magnitude of the velocity is fixed: \v*i\ = v , its direction is updated at the next time step 
by averaging over its neighbors' moving directions: 

1 M 

o t {t + 1) = e(- £fa(t) + + m) (7) 

i=i 

M is the number of neighbors for boid i within radius R: rij — \r\ — fj\ < R. The extra 
interaction term = go(fl — rj)((j^-) 3 — (jr-) 2 )) makes boids repel each other when they 
are closer than l Q , and attract each other otherwise, with l the average distance between 
boids in the flock, this interaction will prevent formation of clusters. The noise term ffi(t) = 
At>(cos(7re;(t)), sin(7re;(t))), where Ci(t) is a random number in the interval [—1,1]. The 
function 0(x) is just the polar angle of the vector x. The position update is simply: r^(i+l) = 
fl(t) + v (cos(9i(t)) , sm(9i(t))) . The parameters in this model are R, 7, Av, v and g . 

The particular form of the interactions should not affect the universal predictions of the 
continuum theory presented above, but rather should only change non-universal phenomeno- 
logical parameters like c, A, D\\ etc. . They also affect the length scale Inl beyond which the 

asymptotic long wavelength forms of the correlation functions (3) and (4) apply. Indeed, an 

si 1 2 

one-loop RG analysis predicts: I^l ~ (IO-D^Dm /AAa) 4=3 x 0(1). Higher loop corrections 
may affect this result, but it presumably remains accurate to factors of 0(1). 

For our numerical model, we estimate (on dimensional grounds): A ~ 1, A ~ (Av ) 2 -^-, 
D\\ ~ D± ~ Inserting these estimates, we find I^l ~ ^-(^^) J=1 - m our simulation, 
choosing units of length and time such that R = t = 1, and taking Av ~ Av c ~ 1/3 in 
these units, for d — 2 we get the lower bound Inl > 30. Previous simulations || took 
Av << 1, and therefore have a much larger I ml- Hence, no non-trivial nonlinear effects 
could be observed since their systems were much smaller than /jvl- 

The above analysis shows that in order to test the scaling behavior with a reasonable 
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system size, one seeks a small Inl by increasing Av and decreasing the radius of interaction 
R as much as possible without entering the disordered phase. In this paper, we report 
the results of a simulation with system size L x L with L = 400 and the number of boids 
iV = 320,000. We choose R — 1, go — 0.6, vo = 1.0, Zo = -707. For these parameter values, 
the flock becomes disordered at Av c ~ .375 as shown in fig. 1(a). The order parameter is 
defined simply as the magnitude of the average velocity of the whole flock:0 = -^| YliLi Vi\- 
To stay in the ordered phase and have enough fluctuations, we choose Av = 0.15. 

Previous simulations have used periodic boundary conditions |§. However, for any finite 
flock, the direction of the average velocity will slowly change, making comparison to the 
analytical results, which assume infinite system size and hence a constant direction for 
< v >, difficult. In order to make < v > constant in its direction, we impose periodic 
boundary conditions in one of the directions, say the x direction, and reflecting boundary 
conditions in the other direction y, i.e., when a boid i with velocity (vf , vf) collides with the 
"walls" at y — ±L/2, its velocity changes to (vf, — vf). The symmetry broken velocity is 
thus forced to lie along the x-direction, without changing the bulk dynamics of the system. 
We will hereafter use "||" and x; "_L" and y interchangeably. 

We first measure the equal time correlation functions. From eq. (3), we predict: 

C p {q) = < 5p{q,t)8p{-q,t) >= J < IM<?»| 2 > ^ 

We see that the equal time correlation function gives us a direct measure of the attenuation. 
The asymptotic behavior of C p (q) can be expressed as: 

C p (q) ~<?r, q±»ql (9) 

~ % K : q± « 4 (io) 

In fig. 1(b), we have plotted the equal time density correlation functions in Fourier 
space: C p (q\\,q± = 2n/L) versus q\\ and C p (q\\ = 0,q±) versus q± from our simulation. The 
scaling behavior at long length scales can be fitted with: C p (q\\,q± = 2n / L) ~ <?|7 2 '° 5 an d 



C P {% = 0,g_i_) ~ qj 1 ' 23 - These two exponents show excellent agreement with the analytical 
results —2 and — | respectively. As can be seen from fig. 1(b), the scaling region for the 
current simulation covers slightly less than one decade in q±. It is not surprising that earlier 
simulations of smaller systems with less carefully chosen parameters (leading to larger Inl), 
did not observe the nontrivial scaling. 

Another interesting measurement of the simulation is the anomalous diffusion of individ- 
ual boids in the direction y perpendicular to the flock's moving direction. We measure the 
"width" of the dispersion of an ensemble of boids: w 2 (t) =< (yi(t) — Ui(0)) 2 >■ The analyti- 
cal behavior of the anomalous diffusion can be obtained from: w 2 (t) ~ J * J * < v y (t')v y (t") > 
dt'dt" where v y (t) is the velocity of the iih boid along y direction at time t. The velocity 
correlation function is given by (4): 

< t>j(0K(f) Vy(X+(f)Xt,t)Vy(X,0) > 

r exp(i(uj - <f>q\\)t)A(u - v s q\\) 2 d 2 qduj _ ^_ 1/c 
J S(q,w) 

which implies: w 2 (t) ~ t 3 ^ 1 ^ = if. In fig. 2(b), we have plotted the width squared w 2 {t) 

versus time i in log- log scale. The scaling can be fitted nicely with w 2 (i) ~ t L3 , which agrees 
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well with the analytical result ts. We have also simulated Vicsek's original model, but with 
parameters Av etc. chosen to make Zjvl as small as possible, and found again w 2 (t) ~ t 13 . 
This supports the universality of our analytic results. 

Besides the scaling behavior, the analytical results (3), (4) also imply the existence of 
sound waves as reflected in the peaks of the correlation functions eqs. (3), (4). From eq. 
(3), at a given value of q, the correlation function has peaks at uj — c±(9 q )q. We have 
measured the power spectrum in the y-direction: < \5p(q\\ = 0, q± = ^-n±,uj = j^n^)] 2 > 
(T=1024) with different values of n±(— 1, 2, • • • , 20). Figure 3 shows the power spectra for 
n± = 5, 10, 20. The spectra are symmetric around u = (we only show half of the spectrum 
for uj > 0) and the positions of the peaks n* versus n y are shown in the inset of figure 3, 
whose slope determines the sound velocity in the y-direction c = 0.62. We have calculated 
the power spectrum of v±_ in the y direction, which shows the same peaks. 
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An interesting phenomenon happens when we calculate the spectrum along the x direc- 
tion, i.e., with q» ^ and q± = 0. As predicted by eqs. (3), (4), we see one single peak for 
each correlation function. Indeed, as shown in figure 4, each power spectrum shows only one 
peak, and again as predicted by (3), (4), the peak for the v± power spectrum is at a different 
uj than the peak of the density power spectrum! This means that the velocity fluctuations 
propagate with a different velocity than the density fluctuations in the x-direction! 

We can then extract from figure 4 the values of v s = 0.93, A = 0.75. (The fact that A ^ 1 
reflects the absence of Galilean invariance). With the value of c = 0.62 determined through 
figure 3, we can predict the sound speeds in all other directions of propagation from eq. (6) 
with no adjustable parameters. To test these predictions, we have also calculated the power 
spectra for the density and the velocity fields at two other angles: tan(# 9 ) = 1/3,4. For the 
large angle 9 q% \ = arctan(4) = 76.0°, the data are shown in fig. 5(a). The peaks for p and 
v± are at the same location, and the wave velocities are c±(8 qj i) = 0.75, —0.37. The data for 
9q,2 — arctan(l/3) = 18.4° are shown in fig. 5(b). The peak at to = c_{6 qt 2) is just barely 
visible in the density correlation, but both peaks show very well in the velocity correlation, 
and the peaks for both correlation functions are at the same locations, giving the velocity 
c ±(#<?,2) — 0.97, 0.59. In fig. 5(c), we have plotted the angle dependence of the wave velocity 
as predicted in eq. (6) in polar angle coordinates (c±(9 q ),9 q ), with the values of v s , A, and 
c determined earlier. We have included in fig. 5(c), the sound velocities for the two angles 
6 q i and 6 q ^- The agreement with the predicted velocities is excellent. 

In summary, the numerical simulations reported here strongly support our analytical 
continuum theory of flocks. The observed sound speeds agree very well with our predictions. 
In particular, our analytical model's assertion that Galilean invariance is absent is confirmed 
by the existence of two different non-zero sound speeds for propagation along the mean 
direction of flock motion. In addition, the sound attenuation shows the anomalous scaling 
we predict 0. 
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FIGURES 

FIG. 1. (a)The order parameter (f>, as defined in the text, versus the noise strength Av. The 
arrow shows the value of Av at which the fluctuations of the ordered state were calculated. (b)The 
scaling behavior of the equal time correlation function for the density fluctuations in the two limits 
as given in eqs. (9,10). The lines illustrate the predicted slopes. 

FIG. 2. The log-log plot of the anomalous transverse diffusion of an individual boid versus time. 

FIG. 3. The power spectrum of the density for different wave vectors. The inset shows the 
peak positions of the power spectrum versus wavenumber. The linear slope determines the sound 
velocity. 

FIG. 4. Power spectra for the density and velocity fluctuations for the same wave vector along 
the parallel direction. The peaks of the two curves are clearly different. 

FIG. 5. The power spectra for the density and the velocity fluctuations in directions 
(a)0 9; i = arctan(4) and (b) 9 q ^ = arctan(l/3). The two peaks are clearly visible, albeit with differ- 
ent magnitudes. In (c), the wave velocities c±(6 q ) are plotted in polar angle coordinates (c±(9 q ),6 q ) 
for the four different directions 9 q = 0, 0^2, tt/2, the two axes represent c x = c±{9 q ) cos (9 q ) 
and c y = c±(9 q ) sm(9 q ) respectively. The solid curve is the prediction from eq. (6) in the text. 
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